Momentum distribution and condensate fraction of a Fermi gas in the BCS-BEC 

crossover 
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By using the diffusion Monte Carlo method we calculate the one- and two-body density matrix 
of an interacting Fermi gas at T = in the BCS-BEC crossover. Results for the momentum 
distribution of the atoms, as obtained from the Fourier transform of the one-body density matrix, 
are reported as a function of the interaction strength. Off-diagonal long-range order in the system 
is investigated through the asymptotic behavior of the two-body density matrix. The condensate 
fraction of fermionic pairs is calculated in the unitary limit and on both sides of the BCS-BEC 
crossover. 

PACS numbers: 



The physics of the crossover from Bardeen-Cooper- 
SchriefFer (BCS) superfluids to molecular Bose-Einstein 
condensates (BEG) in ultracold Fermi gases near a Fes- 
hbach resonance is a very exciting field that has recently 
attracted a lot of interest, both from the experimen- 
tal P, Q and the theoretical side 3]. An important 
experimental achievement is the observation of a con- 
densate of pairs of fermionic atoms on the side of the 
Feshbach resonance where no stable molecules would ex- 
ist in vacuum 0, |^. Although the interpretation of 
the experiment is not straightforward, as it involves an 
out-of-equilibrium projection technique of fermionic pairs 
onto bound molecules 0] , it is believed that these results 
strongly support the existence of a superfluid order pa- 
rameter in the strongly correlated regime on the BGS 
side of the resonance ^ . 

The occurrence of off-diagonal long-range order 
(ODLRO) in interacting systems of bosons and fermions 
was investigated by G.N. Yang in terms of the asymptotic 
behavior of the one- and two-body density matrix [3 . In 
the case of a two-component Fermi gas with spin-up 
and Ni spin-down particles, the one-body density matrix 
(OBDM) for spin-up particles, defined as 



pi(r;,ri) = (v4(r;)V'T(ri)), 



(1) 



does not possess any eigenvalue of order N-^ . This behav- 
ior implies for homogeneous systems the asymptotic con- 
dition pi{r[,ri) ^ as |ri — r'^j 00. In the above ex- 
pression V'|(r) (V'T(r)) denote the creation (annihilation) 
operator of spin-up particles. The same result holds for 
spin-down particles. ODLRO may occur instead in the 
two-body density matrix (TBDM), that is defined as 

P2(r'i,r^,ri,r2) = (V|(rl)V^t(r^)^T(ri)V'i(r2)> ■ (2) 

For an unpolarized gas with N-^ — — N/2, if p2 has an 
eigenvalue of the order of the total number of particles N, 
the TBDM can be written as a spectral decomposition 



separating the largest eigenvalue, 

P2(r'i,r^,ri,r2) = a7V/2^*(r;,r^)^(ri,r2)-f p'2 , (3) 

p'2 containing only eigenvalues of order one. The param- 
eter a < 1 in Eq. ||2J) is interpreted as the condensate 
fraction of pairs, in a similar way as the condensate frac- 
tion of single atoms is derived from the OBDM. 

The spectral decomposition Q yields for homogeneous 
systems the following asymptotic behavior of the TBDM 

P2(r'i,r^,ri,r2) ^ aN/2^*{\r[ - ^D^dn - r2|) , (4) 

if |ri — r'^l, |r2 — rjj 00. The complex function cp 
is proportional to the order parameter {tp^{ri)ipi{r2)) = 
^y aN/2ip{\ri — r2|), whose appearance distinguishes the 
superfluid state of the Fermi gas. Equation should 
be contrasted with the behavior of Bose systems with 
ODLRO, where pi has an eigenvalue of order N Q, and 
consequently the largest eigenvalue of p2 is of the order 
of N^. 

In this Letter we present fixed-node diffusion Monte 
Garlo (FN-DMG) results of pi and p2 for a homoge- 
neous interacting Fermi gas at T = in the BGS-BEG 
crossover. From the Fourier transform of pi{r), we cal- 
culate the momentum distribution of the gas, nk = 
/ (firpi{r)e^^''^ , as a function of the interaction strength. 
From the asymptotic behavior of p2 , we extract the value 
of the condensate fraction of pairs a. The calculated con- 
densate fraction is compared with analytical expansions 
holding on the BEG and BGS side of the Feshbach reso- 
nance. The comparison with mean-field results |^ for 
and a in the crossover region is also discussed. 

We consider a homogeneous two-component unpolar- 
ized Fermi gas described by the Hamiltonian 



H 




(5) 



2 



where m is the mass of the particles and (i', j', ...) 

label spin-up (spin-down) particles. The strength of the 
interaction is assumed to be determined only by the pa- 
rameter 1/kpa, with kp — {Sn'^nY^^ the Fermi wave- 
vector fixed by the atomic density n — N/V, and a the 
s-wave scattering length describing the low-energy colli- 
sions between the two fermionic species. The interatomic 
interactions in Eq. |(SJ are only between atoms with dif- 
ferent spin and are modeled by a short-range potential 
that determines the value of a. In the present study, we 
use an attractive square-well potential, V{r) = —Vq for 
r < Ro and V{r) = otherwise, with hRq = 10~^. We 
have verified that in the density range nRl = 10-'^- 10-5 
the particular form of V{r) is not relevant, and therefore 
the present results are in this sense universal. The differ- 
ent regimes: BEC (a > and l/kpa > 1), BCS (a < 
and l/fci?|a| ^ 1) and unitary limit {l/kpa = 0), are 
obtained by varying the potential depth Vq as in 'lO]. 
Quantum Monte Carlo studies of the Hamiltonian (jSJ 
have already been carried out to investigate the T = 
equation of state 0, ^ and the pairing gap [Tl] . 

In a FN-DMC simulation, the wave function /(R, r) = 
^^(1^)^(1^, t) (R = ri, . . . ,rjv^,ri/, . . . ,rjvj is evolved 
in imaginary time t = it/h according to the time- 
dependent Schrodinger equation, with iprCR-) acting as 
importance sampling function and as nodal constraint. 
The function ^'(R) = ^'(R, r oo), which is the lowest 
energy state of the system having the nodes of V't(R), 
is obtained from the large-time evolution of /(R, r). 
The trial wave function we consider has the general 
form [Ullll ■!/'t(R) = *j(R)*BCs(R), where ^-j con- 
tains Jastrow correlations between all the particles, 

*,j(R) = H /TT(r.,) n fuir^'r)ll fuirw) , (6) 

i<j i^i' 

and the BCS-type wave function ^'bcs is the antisym- 
metrized product of the pair wave functions 4'{vi — r^/) 

*BCS(R) = A {(j){vi - Vv)(j){v2 - V2')..4{'Cn^ ~ r^J) • 

(7) 

The pair orbital 0(r) is chosen as 

0(r)=/3 ^ e^k„.r^^^(^) ^ (8) 

where the sum is performed over the plane wave states 
Icq — 2T:/L{iaxX + d-ayV + ^azz) up to the largest closed 
sheh k^ax = 2n/L{e^^^^ + ^^^^^ + ll.^^^fl'^ occupied 
by iV/2 particles. Here L = T/^/^ is the size of the cubic 
simulation box and I are integer numbers. If 0s (r) = in 
Eq. (jSJ, ^'bcs in Eq- O coincides with the exact wave 
function of a free Fermi gas, i.e., the product of Slater de- 
terminants of spin- up and spin-down particles (l^ . The 
spherically symmetric function 0s (r) in Eq. accounts 
for s-wave pairing. If 1/kpa > —0.2, 0s (r) corresponds 
to the solution of the two-body problem with the po- 
tential V{r), as in Ref. 0. For 1/kpa < —0.2 we use 
instead (j)sir) = 7i[exp(-72r) -1- exp(-72(L - r))]. Here, 
7i: 72, and /3 in Eq. (|SJ| are variational parameters. 



The Jastrow wave function ^j, Eq. ®, is determined 
as follows. For 1/kpa > —0.2, we use f-\i{r) = 1 and 
/tt(^) = flli^) given by the two-body solution of a ficti- 
tious repulsive step potential with range R and scattering 
length a. The boundary conditions /(r — L/2) — 1 and 
f'{r = L/2) = determine the wave function in terms of 
the variational parameters R and a. For 1/kpa < —0.2, 
we use instead f^i{r) = fii{r) = 1 and the model used 
in Ref. |0 for the crossed correlation factor, /|x(r). It 
is worth noticing that the function ipxCR) defined above 
reproduces as a special case the trial wave function used 
in the preceding study |10| , but contains more variational 
parameters. The parameters of the Jastrow function ^'j, 
Eq. ©, are optimized by minimizing the variational ex- 
pectation value {tpT\H\^T) / {iPtH't) ■ The parameters of 
the BCS function ^bcs, Eq. Q, affect the nodal sur- 
face of the trial wave function and they are optimized by 
minimizing the FN-DMC estimate of the energy. For the 
values of 1/kpa used in the present study, the calculated 
FN-DMC energies are in agreement with the results ob- 
tained in although the optimized variational energy 
has significantly improved. 

A direct estimate of any operator O in DMC is known 
as mixed estimate, (0)m = (^t|0|4')/(?At|^), and is ex- 
act only for the Hamiltonian and operators commuting 
with it. If O is a local operator, one can circumvent this 
problem by introducing pure estimators. That is not the 
case for pi and p2 which are the objectives of the present 
work. In order to reduce, and even eliminate in prac- 
tice, a possible bias in the calculation we have used the 
extrapolated estimator (^'|0|^')/(*|*) ~ 2(0)„i - (0)v, 
with (C>)v = (V't|O|V't)/(V't|'0t)- Any residual bias in 
the extrapolated estimator is reduced if the trial function 
VtIR) bas a large overlap with VE'(R). Consequently, 
compared to the calculation of the eigenenergies, the op- 
timization of ipT is a more important issue in the calcu- 
lation of observables hke the OBDM and TBDM [l3|. 

We consider a system with = 66 particles and pe- 
riodic boundary conditions. In Fig. ^ we show results of 
the OBDM pi{r), Eq. for different values of the in- 
teraction strength 1/kpa. In the deep molecular regime, 
1/kpa ^ 1, the OBDM is determined by the molec- 
ular wave function (j)hs{r)- Pi{r) — n/2 J d^r'0Jjg(|r -|- 
r' I )0bs (''')• Eor a zero-range potential the molecular wave 
function is given by 0bs('') oc e^^^"'/r and one finds 
pi{r) ~ ne~'^l°-/2. This behavior is shown in Fig. ^for 
1/kpa — 4. If one moves closer to the resonance, the 
OBDM decays slowly and oscillations start to appear. 
Finally, on the BCS side of the resonance, the OBDM 
becomes more and more similar to the ideal gas result 
pi(r) = ■in{w\[kpr)/{kpr) - cos{kpr)\/{2(kprf\. The 
momentum distribution nt, obtained from the Fourier 
transform of P\(t), is shown in Fig. |21 In the inset of 
Fig. 121 we compare nt, calculated using FN-DMC for 
1/kpa — 4, with the momentum distribution of the atoms 
in the molecular state rik = 4(fcFa)^/[37r(l -I- fc^a^)^] fTsj . 
To reduce finite-size effects in the calculation of the 
Fourier transform for 1/kpa — and —1, we have used 
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FIG. 1: OBDM for different values of the interaction strentgth 
l/kpa (solid lines). The dotted line corresponds to e"'^''" 
for l/kpa — 4 and the dashed line is the OBDM of a non 
interacting gas. 





0.007 
0.006 
1 0.005 












\ 0.004 
\\ 0.003 
;1 0.002 








1 0.001 
1 0.000 








) 1 2 3 4 5 6 









0.0 0.5 1.0 1.5 2.0 

k/kp 



FIG. 2: Momentum distribution rik for different values of 
l/kpa (solid lines). The dashed lines correspond to Wk cal- 
culated using the BCS mean-field approach |Ti| . Inset: nk 
for l/kpa — 4. The dotted line corresponds to the momen- 
tum distribution of the molecular state (see text). It almost 
coincides with the mean-field result (dashed line). 



FIG. 3: Projected TBDM, pf(r), (solid symbols), 
N/2{2pi{r)/nf (lines) and h{r) = pf(r) - N/2{2pi{r)/nf 
(open symbols) for diff'erent values of l/kpa: l/kpa = 1 (solid 
line and triangles), l/kpa = (dashed line and circles) and 
l/kpa — —1 (dotted line and diamonds). 
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FIG. 4: Condensate fraction of pairs a as a function of the 
interaction strength: FN-DMC results (symbols) , Bogoliubov 
quantum depletion of a Bose gas with am ~ 0.6a (dashed 
line), BCS theory using Eq. IIH (dot-dashed line) and self- 
consistent mean-field theory (solid line). 



the following model for the /c-dependence of 



^ 1 



{k/kpf-fi 



(9) 



The values of /.t, A and A are free parameters deter- 
mined by the best fit of the inverse Fourier transform 
of Eq. to the calculated pi(r), with the constraint 
^/Vj2k^^ — For ^ = 1/2, the above expression 

reproduces the standard of BCS theory with /i and 
A, respectively, chemical potential and gap in units of 
the Fermi energy ep = h'^kjp/2m. For l/kpa = and 
-1 the Fourier transform of the BCS model, Eq. (Q, re- 
produces quite well the calculated pi{r) with a /v of 



the order of one. In Fig.|21we also show the results of ?ik 
obtained using the BCS mean-field theory [9|, where the 
values of chemical potential and gap are calculated self- 
consistently through the gap and number equations ^1^. 
As evident from Fig. |21 the mean-field theory sistemati- 
cally overestimates the broadening of the momentum dis- 
tribution, apart from the deep BEC regime where both 
mean-field and FN-DMC results coincide with the mo- 
mentum distribution of the molecular state. 

The condensate fraction of pairs has been obtained 
from the the projected TBDM, defined as [TeLITTt 



P2{r)^j^ / d^rici-'r2P2(ri +r,r2 -f r,ri,r2) . (10) 
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Assuming the eigenvector (p(r) of Eq. (0J normalized 
to 1/V, the large r behavior of P2{r) tends to the con- 
densate fraction a: Vanr^oo P2 {''") — Q^- terms of the 
order parameter, i^(|ri — r2|) — (?/'t(ri)V'i(r2)), one has 
instead \mir^^p^{r) = 2/n J d3r'|F(r')|2. In Fig. 01 we 
show the results of P2{'r)- Finite size effects can be sub- 
stantially reduced if one considers the following decompo- 
sition: limr^ooP2(ri+r, r2 + r, ri,r2) = |F(|ri -r2|)P + 
Pi{r), accounting for the large r behavior of p2 when 
pi(r) is small but not zero. From this result one finds 
for the asymptotic behavior of the projected TBDM: 
limr_oopf (j') = a + N/2{2pi(r)/n)'^. In Fig. El we show 
results for the quantity h{r) ~ pl'(r) — N/2{2pi{r)/n)^. 
For values of l/kpa on the BEC side of the Feshbach 
resonance and up to the unitary limit 1/kpa = 0, the 
asymptotic values of P2(?') and h(r) coincide. On the 
BCS side, instead, finite-size effects become visible in the 
large r behavior of (r) , but are strongly suppressed for 
h(r). The results for the condensate fraction a, as ob- 
tained from the asymptotic behavior of h{r), are shown in 
Fig. El In the BEC regime, the results reproduce the Bo- 
goliubov quantum depletion of a gas of composite bosons 
a = 1 — S-^nmO^/S-^/Tr, where nm = n/2 is the density 
of molecules and = 0.6a is the dimer-dimer scattering 
length jTSIII. In the opposite BCS regime, the conden- 
sate fraction a can be calculated from the result of the 
BCS order parameter holding for r a [Toll 



^BCS('') 



Akp sin(fc^-r) 



(11) 



where ^0 = ^ kpjrat^ is the coherence length and 
Kq{x) is the modified Bessel function. If we include 
the Gorkov-Melik Barkhudarov correction for the pair- 
ing gap A = (2/e)^/3g^g-7r/2fci.|a|^ obtain for 



a = 2/n J c?'^ri^|Qg(r) the dot-dashed line shown in 
Fig. El For Xjkpa < —1 the coherence length ^0 be- 
comes increasingly larger and finite-size effects start to 
be relevant in the FN-DMC estimate of a. For exam- 
ple at 1/kpa = —1, the value of a, as obtained from 
FBCsir) [Eq. (CU], reduces from 0.10 to 0.08 by cut- 
ting off the spacial integral at the simulation box size 
r = L/2. The condensate fraction a can also be esti- 
mated in the crossover region using the self-consistent 
mean-field approach '3| and the result for the order pa- 
rameter: F{r) — l/V'X^k^fe^fc^''""' where Uk and Vk 
are the usual quasiparticlc amplitudes of BCS theory. 
The result is shown in Fig. El with a solid line (see also 
Ref. 0). The mean- field result reproduces the qual- 
itative behavior of the condensate fraction across the 
resonance. However, it does not reproduce the quan- 
tum depletion in the BEC regime, nor the Gorkov-Melik 
Barkhudarov correction to the gap in the BCS regime. 

In conclusion, we have investigated using quantum 
Monte Carlo techniques the one- and two-body density 
matrix of a homogeneous Fermi gas in the BCS-BEC 
crossover. Results for the momentum distribution are 
obtained as a function of the interaction strength. The 
condensate fraction of pairs is calculated from the asymp- 
totic behavior of the two-body density matrix. The re- 
sults obtained are in good agreement with the quantum 
depletion of a weakly-interacting Bosc gas in the BEC 
regime and with the normalization of the BCS order pa- 
rameter in the BCS regime. 
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